Gravitational Ionization: A Chaotic Net in the Kepler System 

C. Chicone* , B. Mashhoon^ and D. G. RetzlofF* 
(February 7, 2008) 

The long term nonlinear dynamics of a Keplerian binary system under the combined influ- 
ences of gravitational radiation damping and external tidal perturbations is analyzed. Gravitational 
radiation reaction leads the binary system towards eventual collapse, while the external periodic per- 
turbations could lead to the ionization of the system via Arnold diffusion. When these two opposing 
tendencies nearly balance each other, interesting chaotic behavior occurs that is briefly studied in 
this paper. It is possible to show that periodic orbits can exist in this system for sufficiently small 
damping. Moreover, we employ the method of averaging to investigate the phenomenon of capture 
into resonance. 



1. INTRODUCTION 

The Kepler problem traditionally involves two point masses mi and m2 moving under the influence of their mutual 
gravitational interaction. The incorporation of relativistic gravitational effects in the Kepler system brings about 
three main post-Newtonian modifications in the traditional picture: 

(i) There exist gravitoelectromagnetic post- Newtonian efl^ects that need to be taken into account. The most impor- 
tant of these is the periastron precession that has played an important role in the historical development of general 
relativity. 

(ii) Moreover, the system emits gravitational radiation and hence there exist radiative post-Newtonian effects in 
the orbit that are characteristic of gravitational radiation damping. This radiative damping is consistent with the 
observed rate of inward spiraling in the Hulse- Taylor binary pulsar system . 

(iii) The system is expected to be affected by gravitational waves that have been emitted by other systems as 
well as by their post-Newtonian tidal perturbations. Gravitational radiation has not yet been directly observed; 
however, the Hulse- Taylor binary pulsar data appear to be consistent with the notion that gravitational radiation 
is emitted by binary systems Half of all stars are expected to be members of binary or multiple systems that 
could emit gravitational waves; therefore, there might be a cosmic background of gravitational radiation that has 
been generated by various sources throughout the history of the universe. Moreover, there could also be primordial 
gravitational waves left over from the epoch at which the Hubble expansion began. The inclusion of all of these effects 
in the Kepler problem is impractical; in fact, the two-body problem in general relativity is intractable. To render 
the problem amenable to mathematical analysis, it is therefore necessary to replace the actual problem by a model 
that contains the main physical effects of interest. We consider a Kepler system in which the main effects of emission 
and absorption of gravitational radiation are taken into consideration. Gravitation is a spin-2 held; therefore, the 
main radiative effects first occur in general at the quadrupole level in emission as well as in absorption. At this level, 
the wavelength of the radiation is large compared to the size of the system. In this work, we analyze the radiative 
perturbations of a Keplerian binary system in the quadrupole approximation. 

In our previous papers [^j4j, we considered the problem of ionization of a Keplerian binary system by a normally 
incident periodic gravitational wave. To render the absorption problem tractable, we simply ignored the emission of 
gravitational radiation by the binary system. Therefore, the influence of the gravitational radiation reaction on the 
binary orbit was not taken into account. In the absence of this dissipative effect, our model of gravitational ionization 
turned out to be a Hamiltonian system to which the basic results of the Kolmogorov-Arnold-Moser (KAM) theory 
could be applied under certain circumstances. The main purpose of the present work is to take gravitational radiation 
damping into account. Thus, we study in this paper the long term nonlinear evolution of a Keplerian binary system 
under the combined action of perturbations due to the emission and absorption of gravitational radiation. 

The nonlinear dynamics of a Keplerian binary that emits and absorbs gravitational radiation is given in our model 

by 

— + — +n^^-elC.,it)r^^ (1) 
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where r := Xi — X2 denotes the relative two-body orbit, r is the length of r, and k 
reaction term TZ is given by 
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where v is the relative velocity, Go is Newton's constant, c is the speed of light in vacuum, and an overdot represents 
differentiation with respect to time, i.e. r = dr/dt. We emphasize here that equation represents the simplest 
equation for the relative motion of two point masses mi and m2 that contains the main dynamical effects of emission 
and absorption of gravitational radiation that we wish to study here; in fact, other post-Newtonian contributions have 
been simply ignored. The justification for this approach emerges from our analysis; that is, the substitution of the 
actual intractable problem by our simple model permits us to arrive at interesting results of possible astrophysical 
interest. 

Let the background spacetime metric be given by g^^ — rj^^ -f- ex^J,v^ where 7/^^ denotes the Minkowski metric, 
e is a small parameter indicative of the strength of the external radiation field (0 < e << 1), and Xtiv represents 
the gravitational radiation field pervading the space. We impose the transverse-traceless gauge condition such that 
Xo^i = 0) (^iXij = 0) s-iid {Xij) is traceless; moreover, Xij satisfies the wave equation ^Xij = 0- Imagine now that the 
mutual gravitational interaction between the two masses is turned off and they are therefore following geodesies of the 
background spacetime manifold. The spatial separation between the two masses is assumed to be small compared to 
the wavelengths of the background (external) waves; therefore, the relative motion of the two masses can be expressed 
via the Jacobi equation. It is useful to express the deviation equation with respect to a Fermi coordinate system 
established along the worldlinc of the center of mass of the system. The Fermi system is the most natural extension 
of a local inertial frame along the path of an observer in spacetime. Hence, the relative acceleration of the two bodies 
in the Fermi system takes on a Newtonian form that is derivable from a quadrupole potential given by ^elCij{t)x^x^ 
with 



K.ij{t) 



^ Xtj 
2 dt^ 



(3) 



where (mi -|- TO2)xcm = miXi -|- m2X2. 

In equation (nf), we have neglected terms of order e^, eS, S"^, or higher, where S is the strength of the radiation 
reaction term as defined in Appendix ^; that is, in conformity with the Newtonian appearance of equation (^ the 
various forces have been linearly superposed. A detailed treatment of equation (jl|) with 7?. = is contained in our 
previous papers , which considered only normally incident external waves for the sake of simplicity. We shall also 
choose here an external periodic monochromatic gravitational wave that is perpendicularly incident on the orbital 
plane. Let this plane be characterized by the vanishing of the third component of r. Then, the nonzero elements of 
the tidal matrix (e/Cjj) are given by 



/Cii = -/C22 = afl'^ cos(fii), 
/C12 = /C21 = cos{ilt + p), 



(4) 



where a and P are of the order of unity and represent the constant amplitudes of the two independent linearly polarized 
components of the incident monochromatic wave of frequency fl and p is the constant phase difference between the 
two components (cf. Q). 

It is interesting to note the partial reciprocity between the emission and absorption of gravitational radiation. In the 
quadrupole approximation, an elliptical Keplerian binary system emits gravitational waves only at frequencies muj, 
m = 1, 2, 3, • • •, where w is the Keplerian frequency of the orbit Q. It follows from a linear perturbation analysis |^ 
that when a Newtonian binary system absorbs gravitational energy from an incident wave of frequency 51, resonances 
occur at 51 = muj, m — 1, 2, • • •. However, it is important to point out that while the binary monotonically loses energy 
in the form of gravitational radiation, the absorption of gravitational radiation energy is not monotonic. The system 
can gain or lose energy in absorption. In our previous papers j^,^, we examined the long term nonlinear evolution of 
the dynamical system with 7?. = and proved the possibility of existence of periodic orbits for which the net flow of 
energy between the binary and the external periodic gravitational radiation field must vanish. In the present work, 
we extend our previous results by taking due account of gravitational radiation damping; in fact, periodic orbits are 
shown to exist for sufficiently small 5/e 

The issue of gravitational ionization provided the original motivation for our work. The possibility of ionization of a 
Keplerian binary by incident gravitational waves had been first discussed within the framework of linear perturbation 
analysis Q. However, the linear perturbation treatment breaks down over time as a consequence of the appearance of 
secular terms in the analysis. On the other hand, the results of the linear analysis can be applied to the possibility of 
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detection — via the solar system — of gravitational pulses [[7| or a stochastic cosmological background of gravitational 
waves ^ ; a review of this general topic can be found in j|] and further work in this direction is contained in jl^-Q and 
the references cited therein. Furthermore, the tidal influence of very long wavelength gravitational waves on galaxies 
has been studied in [^ . Our nonlinear perturbation treatment has led to the conjecture that gravitational 
ionization is tantamount to Arnold diffusion. Moreover, there are interesting circumstances — other than periodic 
orbits — for which ionization does not occur. For instance, it follows from the KAM theorem that when the incident 
wave is circularly polarized (i.e. a = f3 and p = ±7r/2) the perturbed binary is forever confined and ionization does 
not occur if e is below a certain limit (e < Ekam)- It is the purpose of this paper to determine how all of these results 
are affected by the presence of gravitational radiation damping. While our previous work |^,^ has involved time- 
dependent (i.e. nonconservative) Hamiltonian systems with chaotic behavior, the introduction of radiation damping 
leads to a dissipative dynamical system that is not Hamiltonian. 

At each instant of time t, the state of relative motion is determined by the relative position and velocity (r,v); 
however, it is useful to employ instead the six orbital parameters of the osculating ellipse. That is, the unperturbed 
bounded motion is in general a Keplerian ellipse; therefore, it is interesting to describe the state of relative motion 
at each instant by the ellipse that the system would follow if the perturbations were turned off at that instant. The 
orbital elements of the osculating ellipse are closely related to the Delaunay elements employed in this paper. They 
are particularly useful in a Hamiltonian system , since the transformation to the Delaunay action-angle variables 
is canonical and the corresponding generating function is time-independent so that the magnitude of the Hamiltonian 
is unchanged under the transformation. 

A limiting case of the relative osculating ellipse is a circle. In this case, the system emits gravitational waves with 
a single frequency 2w. On the other hand, linear perturbation theory reveals that resonant absorption occurs at wave 
frequencies lu, 2uj, and 3uj. However, it should be remarked that after the passage of time the orbit will no longer 
be circular. The nonlinear perturbation theory employed in our work is based on the ideas and methods originally 
developed by Poincare An important aspect of this theory is its lack of emphasis on the precise form of the 

initial conditions; clearly, this characteristic of the theory is rather useful in celestial mechanics. In this approach, the 
dynamical behavior of a physical system is investigated after transients have died away and a steady state has been 
established. 

In connection with the issue of initial conditions, it is necessary to remark that the external periodic perturbation 
in our model is meant to represent the dominant component of an initial wave packet composed mainly of long- 
wavelength Fourier components consistent with the quadrupole approximation under consideration here. On the 
other hand, high-frequency waves with wavelengths much smaller that the semimajor axis of the system are expected 
to have negligible influence on the dynamics. That is, along any significant portion of the relative orbit the waves are 
expected to go through many oscillations and their net influence is on average expected to be vanishingly small. In our 
model, the strength of the external perturbation is e; at present, it is expected that e ^ 10"^", though gravitational 
waves have not yet been detected in the laboratory. On the other hand, the strength of radiative damping in our 
model is given by 6 (cf. Appendix for instance, for the binary pulsar PSR 1913-1-16, discovered by Hulsc and 
Taylor S ^ 10~^^, while for the Earth-Sun system S ^ 10~^^. The more relativistic a binary system is, the larger 
the magnitude of (5, which for a physical system must always be less than (20\/2)~^ as discussed in Appendix 

Throughout this paper, the external perturbing force is assumed to be due to background gravitational waves. 
However, as explained in detail in our first paper P], this is not necessary since any tidal influence on the binary orbit 
by a far-away mass would contribute a forcing term with the same functional form as in the equation of relative motion 
(|l|). In general, the binary system cannot be isolated from the other masses in the universe; this is simply a result of 
the universal character of gravitational attraction. Thus equation (|l|), which is essentially the damped Newton- Jacobi 
equation, should contain the tidal influence of distant masses as well as background gravitational waves in the forcing 
function. That is, the background spacetime curvature — represented in equation ([^) by the tidal matrix — is due 
to the combined gravitational influence of background waves as well as masses. The analysis of the general problem 
would involve the three-dimensional Kepler system. In our planar model, we neglect the other masses and choose a 
normally incident monochromatic plane wave for the sake of simplicity. 

An interesting feature of equation (|l|) must be noted here: The vector cross product of this equation with the 
relative position vector results in an equation for the rate of change of the orbital angular momentum with time. In 
this result, the torque due to the gravitational radiation reaction is proportional to the orbital angular momentum 
according to equation (^; therefore, the action of this torque is only to change the angular momentum normal to 
the orbital plane. The torque due to the external force is similar in character, since the external wave is transverse 
and is assumed to be normally incident on the orbital plane. The net result is that the relative motion is planar, 
and this fact considerably simplifies the analysis of the relative motion. Moreover, in conformity with the quadrupole 
approximation the net velocity of the center of mass (i.e. the velocity of the binary system as a whole) is uniform. 
Therefore, we choose the "inertial" reference frame (t, x) for our model such that it is comoving with the whole binary 
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system and that Xcm = 0. In this frame, the orientation of the coordinate system is so chosen that the perturbed 
Keplerian motion has a positive initial orbital angular momentum. 

The motion in our model is planar; therefore, it is interesting to express equation (|^) in terms of polar coordinates 
(r, 6) in the orbital plane. Thus, 

dr 

dO^ _ Pe_ 
dt~^' 

dpr ^ k ^ pI ^ l&Glmim2 Pr f 2^gP|^ffc\ 
dt 5c^ \ ^ 3 r y 

-erU^ [a cos 26 cos Vlt ^- (3 sin 26 cos {Vlt + p)] , 
dm^SGlm,rn,pef pJ^^k\ 
dt 5c° \ r J 

+er^9? [a sin 26 cos Vtt - (3 cos 29 cos [Vlt + p)] . (5) 

Here pe is the orbital angular momentum and the initial conditions are chosen such that > 0. It is clear from these 
equations that the corresponding vector field in phase space is periodic with frequency fi. 

The main dynamical result of our work may now be stated qualitatively in terms of the characteristic forms of 
the perturbing functions. The gravitational radiation reaction force is proportional to r~'^ with q > 3. Thus if r 
becomes small, radiative damping takes over and leads the system inexorably to collapse. On the other hand, the 
external potential is proportional to — i.e. the energy exchange is proportional to the orbital area through which 
the normally incident wave passes (cf . Q ) — so that if r becomes large the external force takes over and the results 
of our previous work ||^,^ imply that Arnold diffusion might take place leading to the ionization of the binary system. 
Thus it appears that the long time dynamical behavior of the system under consideration here is characterized by 
two possibilities: collapse to the origin and unbounded growth; in fact, our results indicate that the collapse scenario 
is prevalent in most cases. 

The plan of this paper is as follows: In Section 2, we discuss the radiation reaction term TZ in equation (^. Section 
3 is devoted to the determination of the bifurcation function and the proof of the existence of periodic orbits in the 
damped system. In Section 4, the average behavior of the dynamical system (|l|) with e = is determined for all 
times. Section 5 describes the average properties of the system (^. Section 6 considers the case of circularly polarized 
incident waves. Section 7 describes the chaotic net using numerical analysis. A number of detailed computations are 
relegated to the Appendices. The paper relies on our previous work [^,^ for background material regarding the topic 
of gravitational ionization; however, we have tried to make this paper essentially self-contained. 



2. GRAVITATIONAL RADIATION DAMPING 



In the quadrupole approximation under consideration here, gravitational waves carry energy and angular momentum 
away from the system but not linear momentum. Similarly, in absorption the orbit exchanges energy and angular 
momentum with the incident gravitational wave but not linear momentum. This implies that the motion of the system 
as a whole — i.e. the center-of-mass motion — is not affected by the emission and absorption of radiation in the 
quadrupole approximation. Taking proper account of the emission of gravitational radiation, the conservation laws of 
energy and momentum require that these losses be reflected in the motion of the system. To satisfy this requirement, 
it is sufficient to include in the equations of motion of a particle in the system a gravitational radiation reaction force 
per unit inertial mass given by A, 

where {Vij) is the quadrupole moment of the system (cf. [0 and the references cited therein). We may look upon 
this force as the spin-2 analog of the standard spin-1 radiation reaction force in electrodynamics JlSf . The radiation 
reaction term in the equations of motion cannot be derived from a potential; otherwise, the dissipative dynamical 
system would be Hamiltonian. 

Let £ denote the energy radiated as gravitational waves; then, 

d£ Go ij 
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where 



= J p{3x'x^ - 5.,jX^)dV. 

Here the quadrupolc moment of the system is traceless by definition, p is the mass density, and dV is the volume 
element. Equation (|^) is a standard result of general relativity in the quadrupole approximation; that is, it can 
be obtained from a detailed treatment of the linearized form of the gravitational field equations together with the 
physical interpretation of the Landau-Lifshitz pseudotensor. Following this approach, a similar expression can be 
derived in the quadrupole approximation for the angular momentum carried away by the gravitational waves. The 
energy and angular momentum radiated away via gravitational waves are lost by the orbit; therefore, the equation of 
orbital motion should reflect this loss. For both energy and angular momentum, this must be accomplished by the 
introduction of the radiation reaction force. By analogy with electrodynamics, we expect that the radiation reaction 
force per unit inertial mass is of the form given by equation (|6|) , so that for the binary system under consideration we 
can express the equations of motion as 

d^x\ Gom2(xi -X2)' _ 2Go d'>V,j ^ 
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where r = Xi — X2. Multiplying the first equation by mi and the second equation by TO2 and adding them results in 
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In the absence of perturbations, Xcm is at the origin of the coordinate system; therefore, this equation of motion for 
implies that it is effectively unchanged in our approximate treatment. It is thus possible to set Xcm = 0, so 
that the emission or absorption of gravitational waves does not affect the total linear momentum of the system in the 
quadrupole approximation. The relative motion is obtained by subtracting the two equations: 



V ^ Go(mi+m2)r' _ 2Go d^V,, ^, 



df^ r3 15c5 dt^ 



eJC,,ity. (10) 



It follows that the full content of the equations of motion is simply contained in the expressions xi = m2r/M and 
X2 — — mir/Af , where M = mi + 1112 is the total mass and r is given by the equation of relative motion (RW- The 
quadrupole moment of the system is then given by 

Vij — mi{3x\xi — SijXi) + 7712(8x2X2 — SijX^) — p,{3r''r^ — (11) 

where p, is the reduced mass, i.e. p, — m-\milM . 

It is clear from the form of equations (p^ ) and ( p] ) that the relative motion must be determined from an equation 
in which the order of differentiation with respect to time exceeds two. On the other hand, the higher derivative 
terms only appear in the radiation reaction acceleration that has an overall dimensionless strength given by a small 
parameter b that is particularly small for Keplerian (i.e. nonrelativistic) binaries as discussed in Appendix ^ It 
follows that the problem can be avoided in this case if we proceed iteratively and substitute for from equation ( |l0| ) 
every time it appears in the calculation of gravitational radiation damping. Then, keeping only terms of the first 
order in the perturbation parameters e and 8 in equation ( |l0| ) we find that 

1 d^T>ii kr f ^ n 8fc\ „ kf , ^ r,- - 

^ = -2-r 9v^ - 15r^ % - 90— {3v^ - If^'r^ 



p, dt^ \ r 
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which is the approximate expression to be used in the radiative damping term 
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in order to obtain equation (g). In deriving equation (|12|), we have repeatedly used the fact that r • v = rf. In the 
following sections, we shall use equations (0) and (^) to study the dynamical behavior of the Kepler system over a 
long period of time. 

Let us now return to equation ( |l0|) and point out that this equation of relative motion — in which the effects of 
emission and absorption of gravitational radiation have been taken into account in the quadrupole approximation and 
all other retardation effects have been neglected — is consistent with the conservation laws of energy and momentum. 
We will illustrate this point explicitly for the energy of the system; the case of angular momentum is similar but will 
not be treated here. Let us therefore derive an energy expression for the equation of relative motion (|l^). Multiplying 
both sides of that equation by the reduced mass and the relative velocity results in 



d 
dt 
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45^ 



4:5c 



(14) 



where we have used Vij = /i [3(uV-' + r'^v^) — 2r • v Sij] and the fact that ) is symmetric and traceless. The first 
term in the square brackets in equation ( p^ ) is the Keplerian energy of the relative orbit 2E = iiv^ — 2kfi/r, 
and the second term is a relativistic contribution. The rate at which this latter term varies averages out — over the 
period of the unperturbed Keplerian motion — to a quantity that is negligible at the level of approximation under 
consideration in this paper. It follows that the average rate of loss of Keplerian energy of the dynamical system by 
radiation damping is equal to the average rate at which energy leaves the orbit via gravitational waves, as expected. 
This latter rate is given by averaging equation (|^) over one period of elliptical Keplerian motion and the result is ||^ 
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where a and e are the semimajor axis and the eccentricity of the Keplerian ellipse, respectively. At this average rate, 
gravitational radiation energy permanently leaves the system and goes off to infinity. No energy is actually lost; the 
orbital energy is simply converted into gravitational radiation energy. 

It is sometimes stated (cf. |l^ ] and references therein) that the force due to gravitational radiation damping could 
be derived from a quadrupole potential. If so, the dynamical system under consideration here would be Hamiltonian. 
That is, a simple comparison of the two perturbing influences in equation (|l(j) reveals that this system would be 



Hamiltonian if d^Vij/dt^ were simply a function of time. However, this cannot be the case since d^T>ij/dt^ is a 
function of r and v; moreover, the system is manifestly dissipative. 



3. CONTINUATION OF PERIODIC ORBITS 



The investigation of the dynamics of the perturbed Kepler system (|l|) can be simplified considerably if the equations 
of relative motion are expressed in terms of action-angle variables that are suited to the unperturbed system ||l^,0 . If 
the perturbations of the Kepler system ^ are turned off at a given instant, the relative motion will follow an osculating 
elliptical orbit with semimajor axis a, eccentricity e, eccentric anomaly u, and true anomaly v. The orientation of 
the underlying coordinate system can be so chosen that the relative orbital angular momentum is positive. The 
action-angle variables appropriate to the Kepler problem are the Delaunay elements [^,0 that are given for the 
planar problem under consideration by 
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To analyze the dynamical behavior contained in equation 
equations in dimensionless form. This is done in Appendix 

expressed in terms of Delaunay elements; this transformation is presented in Appendix 
form of the dynamical equations: 
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and fLjfoj fe^ fg can be expressed in terms of Delaunay elements using classical methods of celestial mechanics 
involving the Bessel functions (cf. Appendix |b|) . The right hand side of the system ( p^ ) is periodic in time with 
frequency fl. For 5 = 0, the system ( p^ ) has periodic orbits as proved in moreover, this Hamiltonian system 
exhibits dynamical behavior |^,^ that appears to be characteristic of Arnold diffusion 23|. In the following, we 
show that when radiation reaction is taken into account the periodic orbits persist for sufficiently small radiative 
damping. For e = 0, on the other hand, the system (^) continuously loses energy to gravitational radiation and 
eventual collapse is inevitable; moreover, this dynamical behavior is structurally stable and is therefore expected to 
persist for sufficiently small e ^ 0. 

Let us now assume that there are relatively prime positive integers m and n such that muj = nQ, i.e. the unperturbed 
Keplerian orbit is in resonance with the disturbing function. Then, it turns out that the three-dimensional resonance 
manifold (period manifold) 



is a normally nondegenerate manifold 
is given by 



= {(L, G, £, g) : muj nfl} 
Moreover, the solution of the unperturbed system starting at (Lq, Gq, io, go) 

t i-> {Lo,Go,ujt + £o,go), 



where t = t — to and to is the initial time. It turns out that we can set to = with no loss of generality. It follows 
from 1^,^,^ that we must project the partial derivative of the Poincare map with respect to e onto the complement 
of the range of the infinitesimal displacement. The partial derivative of the Poincare map with respect to e at e = 
is obtained from the solution of the second variational initial value problem 
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S/e so that < A < oo. In fact, the partial derivative is given by the vector 

27r 



(19) 



(20) 



27r 

m—,Lo,GoJo,go 



, Gf 



m—,Lo,Go,£o,go 



7 



4 ^TO^,io,Go,^o,.go^ ,ge (^m^,LQ,Go,£Q,go 

Here the independent variables are (t, L, G, i, g) and the last four variables of each component function give the initial 
point for the original solution of the unperturbed system, i.e. equations ( pT[ ) with e = and 5 = 0. 

The appropriate projection is onto the first, second and fourth coordinates. In fact, the bifurcation function is this 
projection restricted to , and is given by 
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where L = Lq is fixed by the choice of resonance, i.e. m/ Lq = nfl, and Go > by our choice of spatial coordinate 
axes. Dropping subscripts indicating initial values of Delaunay elements for the sake of simplicity, the components of 
the bifurcation function ( pT[ ) can now be written as 
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where X is defined by 



and C and S are given in Appendix The integral in (|2^) has been calculated in with the result that X = for 
n > 1, while for n = 1, 

X = —TTma^illa^Arn cos m£ cos 2g — Bm sin m£ sin 2g) 

+f3[A,n cos {m£ — p) sin 2g + B,n sin {■m£ — p) cos 25]}. (24) 

The periodic orbits that continue are determined by the simple zeros of the bifurcation function ( p^ ) . The integrals 
multiplying A in equations ( |2^ ) are given in Appendix therefore, periodic orbits continue to exist when the following 
equations are satisfied by simple zeros: 

dX ^ ^ /„ 73 2 37 4\ 

^ + 27rA4T(8 + 7e2) = 0, 
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For n > 1, X = and hence these equations have no solution as the expressions multiplying A are manifestly positive. 
For n = 1 and A = 0, we have proved in |^ that there are simple zeros of the bifurcation function for all positive 
integers m. It follows that equation ( |25|) must have simple zeros as well for sufficiently small A. 

It is important to remark here that this conclusion would not be altered by the inclusion of terms of order e^, eS, 
5"^, or higher in the original system (^. In fact, the incorporation of such higher order effects could only affect the 
shape of a periodic orbit of the type investigated here but not its existence. 

The Kepler system under investigation in this paper may in general have other periodic orbits than those revealed 
by our first order method. The investigation of the stability of the periodic orbits under consideration is beyond the 
scope of this work. If there is a stable periodic orbit for e > and S > 0, then its (open) basin of attraction consists 
of trajectories that are permanently captured into resonance. This same issue will be discussed in connection with 
the averaging methods that we introduce in the next section. 
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4. AVERAGED DYNAMICS 



In this section, we wish to consider our system in Delaunay elements when no external gravitational waves are 
present. The motivation for this discussion is the experimental observation of the inward spiraling of the members of 
the Hulse- Taylor binary pulsar j^] . Thus we consider the traditional Kepler problem except that the post-Newtonian 
effect of gravitational radiation damping is taken into account in the quadrupole approximation. We expect that the 
system remains planar and loses energy and angular momentum such that eventually L — > 0, G ^ 0, and the system 
collapses. The behavior of the system in the infinite past is less trivial. The precise manner in which the orbit behaves 
on average for t — > ±oo can be determined using the method of averaging. In this case, we have 

L = 6fL, G = SfG, 
i^uj + Sh, g = 5fg. (26) 

Equations ( p6|) are of the general form 

/ = e/(/,^,e), ^^Lu{I) + eg{I,v,e), (27) 

where both / and g are 27r periodic in the angle variable ip. In fact, in our case / = (L, G, g) and (p = i. The Averaging 
Theorem Q asserts the following: Suppose that t t-^ (lit), ip{t)) is a solution of equations ( p7| ) and Lu{I(t)) > for 
all t. If 1 1— > J(i) is the solution of the averaged system 

r2Tr 



with J(0) = /(O), then, for sufficiently small e, there is a constant C independent of e such that 

\I{t) - J{t)\ <Ce for 0<t<-. 

e 



That is, the "action variables" of the original system ( |27D remain close to the solution of the averaged system over a 
long time-scale. 

It is important to note that the Averaging Theorem, in the generality stated above, is only true for the case of a 
single frequency, i.e. there is only one angle variable and its frequency does not vanish. If there are more than one 
angle present, then resonances among the frequencies must be taken into account. Though g is an angle variable, its 
evolution is slow and its average rate of variation vanishes; hence, it may be regarded as an "action" variable for the 
purposes of this section. Therefore in our case, the averaged system obtained from (|26|), where we use L, G, and g 
for the averaged variables, e = \J L?- ~ / and IS. — 6 / e, is given by (cf. Appendix O 

■ A / 73 2 37 4 



G7 V 3 12 
A 



G = -e-^(8 + 7e2 



5 = 0. (28) 

To determine the dynamics of the averaged system (|^ for the variables L and G, we note that equations ( |2^ ) 
are autonomous. Of course, the system is in this case singular on the (invariant) L-axis. Also, recall that it suffices 
to consider only the case where G > 0. After multiplication by 12L^G^/eA, we obtain the dynamically equivalent 
system of equations 

L = -LiSlC^ - SeeG^L^ + 425^"*), 

G = -G3(180L2 - 84G2). (29) 

The Delaunay elements are defined only for L > and < G < i. It is easy to determine the phase portrait of the 
system (|29|); in fact, all orbits in the sector bounded by the L-axis and the line L = G are attracted to the origin. The 
"singular set" corresponding to the limits of the definition of Delaunay elements, e = and e = 1, consists of the line 
L ~ G and the L-axis, respectively; indeed, these are invariant sets in the scaled system (kS). Recall that L = ^/a 



and G = ^ a(l — e^), where a and e are the semimajor axis and the eccentricity, respectively. From our analysis, it 
follows that a approaches zero as i ^ oo and a increases without bound as t — > — oo. We claim that e ^ as i — > oo 
and e^last^— oo. To see this just note that 
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= ^(^2 _ g2)(425L2 - 121G2); 
now let r] = G/L and observe that < 77 < 1, 77 = y/1 — e"^, and 

By rescaling time, the dynamics of this equation turns out to be equivalent to the dynamics of 

T]' = r]{l - ?7)(425 - 121r/2), 

where a prime denotes differentiation with respect to the new temporal variable. The last differential equation has 
a source at 77 = 0, a sink at 77 = 1, and no other rest point on the interval (0, 1) on the corresponding phase line; 
therefore, we see that 77 — > 1 as < ^ 00 and 77 ^ as t ^ —00, i.e. e ^ as t 00 and e ^ 1 as i — > —00. Hence 
the orbit is unbound in the infinite past. 

The main physical conclusion of this section has been previously obtained by Walker and Will [0. However, the 
method of averaging employed here provides a simple and transparent proof of their results on the dissipative nature 
of the radiation damping. 



5. PASSAGE THROUGH RESONANCE 

In this section, we consider some additional aspects of the dynamics of equation (|^) that can be determined using 
the method of averaging. In particular, we will discuss the phenomenon of capture into resonance. 

The Delaunay "act ion- angle" variables provide the appropriate form of the dynamical equations required for the 
averaging method. Here we will begin with system ( ]l7| ) in the form 

L 
G 
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dH* 




^ de 


dn* 


+ eA/G 


dg 


dH* 


+ eA/g, 


' dG 


1 


dn* 






n, 





^ = — + + eA/£, 

(30) 
where 

n* ^ ^r^^ [aC{L, G, e, g) cos s + 13S{L, G, £, g) cos(s + p)] , 

denotes the 0(e) terms of the Hamiltonian giving the gravitational wave interaction in Delaunay elements and s := Qt. 
Here we assume that e is the small parameter; however, our conclusions remain the same if 6 — eA is taken to be the 
small parameter in this system. 

A key observation that motivates the analysis of this section is the fact that the angular variable g is a slow variable 
for (|30|). Thus, the dynamical system can be treated as a two-frequency system with fast variables £ and s. The 
associated resonances are given by relations between the frequencies and f2. More precisely, if m and n are 

relatively prime integers, the (m : n) resonant manifold is given by 

{{L,G,gJ,s) : = nVl}. 

After averaging system (|3^) over the fast variables, we obtain the system (|2^) that was analyzed in the last section. 
Thus, according to the Averaging Principle [^, every perturbed Keplerian elliptical orbit collapses even in the 
presence of gravitational wave interaction. However, as is well known, the Averaging Principle is not valid for two- 
frequency systems unless additional hypotheses are imposed. Violations of the asymptotic estimates embodied in the 
Averaging Principle are associated with the dynamical phenomenon called capture into resonance. A theorem of A. 
I. Neishtadt (cf. |Q, p. 163, for a discussion) allows for the possibility of capture into resonance and provides a strong 
result on the applicability of the Averaging Principle. 
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Neishtadt's Theorem requires two hypotheses that we refer to as condition A'' and condition B. Condition N requires 
that the rate of change of the frequency ratio of the fast angles with respect to time along the averaged system be 
bounded away from zero. For our case, the frequency ratio is given by (l/L'^)/f2, and the required derivative is 



Clearly, the derivative is bounded away from zero as long as L and G are bounded away from infinity. As we have 
shown in Section ^, both L and G decrease with time in the averaged motion. Thus, condition N is satisfied along 
each fixed orbit of the averaged system. Condition B is generically satisfied for most orbits. 

If conditions N and B both hold, then Neishtadt's Theorem asserts that for all initial points outside a set of measure 
not exceeding a constant multiple of ^/e, the averaged motion approximates the motion over a time-scale 1/e with 
an error given by 0(| lne|-\/e). In particular, capture into resonance is rare, and most motions are well approximated 
by the averaged motion. This is the behavior predicted for our model equations. That is, except for a set of initial 
conditions with small measure, the variables L and G associated with the perturbed orbit are such that L — s- and 
G ^ as < ^ oo; the long time behavior of the osculating ellipse is described by eventual collapse. 

The exceptional set of initial conditions mentioned in Neishtadt's Theorem is not empty in our case. In fact, we 
have proved that some of the unperturbed periodic orbits can be continued to periodic orbits in the presence of 
perturbation for sufficiently small e and A; indeed, these periodic orbits are in a sense permanently captured into 
resonance. 

In order to determine the average dynamics in more detail, the dynamical behavior near the resonant manifolds 
must be studied. One approach to the study of this behavior is provided by the theory of partial averaging near 
a resonant manifold. In fact, these ideas are used in the proof of Neishtadt's Theorem. An elementary exposition 
of the main aspects of the now well established theory can be found, for example, in p^ ; a mathematical proof 
of Neishtadt's Theorem is given in psf . This theory requires second order averaging near the resonant manifold, 
where we find essentially a pendulum-like equation with slowly varying amplitude, phase, damping, and torque. The 
existence and stability of its stationary solutions together with the positions of the associated invariant manifolds 
determine the actual average dynamics near each resonance. We have obtained the second order partially averaged 
system for the perturbed Kepler problem, but its analysis is beyond the scope of this paper. However, we mention 
here two facts in this connection: Near (m : n) resonance with n 7^ 1, there is no capture into resonance; moreover, 
as A increases, the likelihood of capture into resonance decreases. 

Figure |l| shows an example of capture into resonance. The top panel in this figure illustrates capture into resonance 
and the middle panel involves an orbit whose initial conditions are slightly perturbed compared to the orbit depicted 
in the top panel. The middle panel shows passage through resonance on the time-scale of our numerical experiment. In 
the bottom panel, the behavior of the orbital angular momentum G is plotted versus L = yfa for the system depicted 
in the middle panel. During this passage through (1 : 1) resonance, the energy of the orbit is constant on average 
while the orbital angular momentum increases in this case — i.e. the eccentricity of the osculating ellipse decreases 
from e k, 0.8 to e ~ 0.4. That is, there is balance on the average between the processes of emission and absorption 
for energy but not for angular momentum. Recall that the resonance condition only fixes the energy of the osculating 
orbit and so the angular momentum can change. A preliminary analysis indicates that this behavior is consistent with 
the second order partially averaged dynamics; however, a full investigation of this interesting phenomenon necessitates 
further research. It would be quite interesting, of course, if this theoretically rare phenomenon could be observed 
astronomically: A binary system gradually spirals inward, but when the decreasing semimajor axis reaches a certain 
value corresponding to resonance with an external periodic perturbation the collapse process temporarily ceases — 
though the binary orbit's eccentricity could change considerably in this sojourn while the semimajor axis fluctuates 
with increasing amplitude about the resonance value — during the period of passage through resonance that may be 
very long compared to the orbital period until the collapse process resumes again. We also note that for our choice 
of e the time-scale for the validity of the averaged motion for orbits not captured into resonance is 1/e = 10*, a value 
that is an order of magnitude smaller than the integration time of 10^. Finally, we emphasize that the numerical 
experiments are conducted by integrating the original equations of motion — not the averaged system. 



A remarkable outcome of our nonlinear analysis of gravitational ionization has been the recognition that for 
normally incident monochromatic plane waves with definite helicity (i.e. right or left circularly polarized waves) the 
KAM theorem would ensure that there would be no ionization for e < Ckam- This result has been further discussed 
in B. It is clear on physical grounds that the inclusion of radiative dissipation cannot change the basic result of 




(31) 



6. CIRCULARLY POLARIZED FORCING FUNCTION 
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this analysis: The system would constantly lose energy and final collapse would be inevitable. While the KAM 
theorem only guarantees confinement, damping by radiation cannot but lead to catastrophic collapse. For e > Ckam, 
the situation is less clear and — as discussed in previous sections — there could be an interesting interplay between 
collapse and Arnold diffusion. 

It must be pointed out that these results are independent of the frequency of incident radiation and only require 
that a ~ f3 and p = ±7r/2. The dynamical system can be viewed from a frame rotating with frequency 51/2; in this 
frame, there exists an energy integral once the radiative damping is neglected. That is, the incident gravitational wave 
stands completely still as a consequence of the gravitational helicity-rotation coupling [26| and hence the dynamical 
system is time independent in the absence of damping. The KAM theorem would then imply that for sufficiently small 
e the system is forever bounded. This confinement brings into rather sharp focus the nonreciprocal nature of emission 
and absorption of gravitational radiation. Clearly, the continuous loss of energy by the binary would eventually result 
in the collapse of the system. 



7. TRANSIENT CHAOS 

The equations of motion can be integrated numerically given a set of initial conditions at t — 0. For this 
purpose it is convenient to express these equations in dimensionless form using the scale transformations discussed in 
Appendix The actual system that we use is given by 

d9^ _ P9_ 

-erfi^ [a cos 20 cos ilt + /3 sin 29 cos (fit + p)] , 

+er'^Q'^ [a sin 26 cos Vlt - f] cos 20 cos {Qt + p)] . (32) 

We remark that this system is numerically ill-conditioned near collision, i.e. near r — 0. Unfortunately, there does 
not seem to be a way to regularize the system for numerical integration as can be done for regular perturbations of 
Kepler motion as discussed in p7[ | . The problem is that the two perturbation terms in our system enter with negative 
and positive powers of the variable r. 

A possible set of initial conditions for ( p2[ ) that we use for our numerical experiments in this work is {pr,pe, r, 6) — 
(e, 1, 1,0), which corresponds to an osculating ellipse at t = with unit orbital angular momentum, eccentricity e, 
semimajor axis a = (1 — e^)~^, and Keplerian frequency w = (1 — e^)^/^. The true anomaly of this osculating ellipse 
is given by w = 7r/2, the Delaunay element g — — 7r/2, and the periastron occurs aX 9 = —ix 12. 

The main conclusion of the analysis in this paper is that for sufficiently small e, the majority of initial conditions 
lead on average to collapse. However, from basic results in nonlinear dynamics, we expect that near resonances there 
would be transverse intersections of stable and unstable manifolds of some of the perturbed periodic orbits. This 
suggests the existence of transient chaos — chaotic invariant sets that are not attr actors. We will see in a moment that 
numerical experiments are consistent with this observation. However, since we do not have estimates on the size of the 
perturbation parameters for which the expected behavior is valid, there does not seem to be a mathematical obstruction 
to the existence of attracting chaotic sets (strange attractors) for some choices of the perturbation parameters. In 
this regard, we mention the connection between the planar Kepler system and a two-dimensional anharmonic system 
with an external periodic force and damping ||2^. The existence of transient chaos for small amplitude forcing and 
the possibility of the existence of strange attractors for larger amplitude forcing and damping is an interesting and 
well documented dynamical feature of the systems in the latter class. In fact, the existence of transient chaos and 
strange attractors for the case of a single oscillator is well known (see, for example, [^). It should also be mentioned 
that the connection between the Kepler system and oscillators is quite general; in this paper, we have restricted our 
attention to the planar Kepler problem for the sake of simplicity. It is possible that the three-dimensional Kepler 
system given by the general form of equations (|^) and (^ would exhibit novel features not found in the planar case. 
We have conducted a 'preliminary numerical search for a strange attractor in our perturbed Kepler system, but we 
have not succeeded in finding such an attractor. This is an interesting problem for further research. 

As a numerical example consistent with transient chaotic behavior, we refer to Figure ^ where the evolution of 
the action variables, strobed at each cycle of the incident wave, is depicted. While both action variables tend to 
eventual collapse, complicated transient dynamics is indicated in this case. It is interesting to note that here the orbit 



dr 
'di 

dpr 

~dt 

dpe 
dt 
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passes through a dense set of (m : n) resonances that may be responsible for the comphcated behavior of these action 
variables. The evolution of these variables in the bottom panel should be compared and contrasted with that in the 
bottom panel of Figure 0. 

Finally, we remark that while energy continuously leaves the system via the emission of gravitational waves, the 
definite trend toward collapse that results may — under special circumstances — be counter-balanced by the input 
of energy into the system by the external periodic perturbations. If the net flow of energy into the system is positive 
on average, gravitational ionization will take place — the semimajor axis of the osculating ellipse will grow, albeit 
not necessarily monotonically. This balance seems to depend very sensitively on the choice of parameters and initial 
conditions. For instance, let the initial conditions be {pr,pe,r,6) = (0.5,1,1,0) in the system ( p2[ ) and set the 
parameters such that a = 2, /3 = 2, p = 0, and 17 = 1.299038106; for e = 0.005 and 5/e = 0.00198992 the system 
appears to ionize over a long time-scale, while for 5/e — 0.00198993 the system collapses. 



APPENDIX A: SCALE TRANSFORMATIONS 



To facilitate the analysis of the dynamical behavior of the Kepler problem with both incident gravitational waves 
and gravitational radiation damping, we transform the equations of motion for this system to dimensionless form. To 
this end, let = R^f^ and t = T^t, where Rq and Tq are scale parameters and and t are dimensionless. This will 
be true for all quantities with a hat. Substituting for r* and t in equations (0) and (||) leads to 



'2 



I2v^ - 30f 



-50f 




eJCijP = 0, 



(Al) 



where 



k — 



4 G'^mim2 

5 c5Toi?o 



and 



Let us now consider the physical interpretation of such a transformation. If we take a binary system with initial 
period 27rTo and semimajor axis i?o, then by Kepler's law 



Go (mi + m2) 



1 



i.e. k = 1. In general, one can define uj — cjTo and a — a/Ro; then, Kepler's law becomes cD^ = k/d^, where a 
is the dimensionless semimajor axis. We also note that 5 is dimensionless and turns out to be a small parameter, 
< 5 « 1, for all realistic binary systems. By setting A; = 1 and dropping all hats in equation dAlt), we obtain 



12v^ - 30r^ 



r I dt r 



-50f^ + — 
3r 



eJCijT^ = 0, 



(A2) 



which is the form of the equation of motion that we investigate in this paper. Now everything is dimensionless in 



equation (A_2); there are two small parameters e and 5 and because k — 1, the underlying scales i?o and Tq are 



connected such that i?Q — kT^ . Furthermore, if we take the direction of incidence of the external gravitational waves 



to be perpendicular to the initial orbital plane, then it follows directly from (A2) that the damped motion is planar. 

Let us now estimate the magnitude of 5 for binaries of physical interest. Expressing 5 in terms of the total mass 
M = TOi -|- TO2 and the reduced mass /i = TOim2/M, we find 



<5 = lfii 
5 \M 



Ro \ ( GqM 



Furthermore, our assumption that fc = 1 implies 
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cTq \ c^Ro 

hence, we have 

1 fGpMY'' 
5{mJ{c^RoJ " 

Here Aji/M < 1 by definition, and GqM/c^Rq is always less than 1/2. The quantity 2GoM / c^Rq is the ratio of the 
Schwarzschild radius of the system to the semimajor axis of the binary; thus, the system is nearly a black hole for 
2GqM I c^Rq ~ 1, while for a Keplerian binary 2GqM /c^Rq « 1, which indicates that the semimajor axis of the orbit 
is much larger than the gravitational radius of the binary system. It follows from these remarks that for a physical 
system 5 < 0.04. For the Earth-Sun system, for instance, fi/M ~ 3 x 10"*^ and GoM/c^Rq ~ 10"^, so that S ~ lO"^'^; 
by comparison, S is negligibly small for an artificial satellite in orbit about the Earth. On the other hand, for the 
Hulse-Taylor binary pulsar, PSR 1913+16, we find that S ~ lO'^^ 




APPENDIX B: EQUATIONS OF MOTION IN DELAUNAY 

ELEMENTS 



Let us now transform the equations of motion (|A2| ) from Cartesian coordinates and velocities to Delaunay elements. 
The standard derivation of the equations of motion in Delaunay variables is based on the assumption that the system 
under consideration is Hamiltonian. However, the dynamical system in our treatment is dissipative; therefore, in this 
appendix we present a ge neral derivation for the planar Kepler problem. This transformation is performed in two 
steps. First equation ([A2|) is expressed in polar coordinates and then in Delaunay elements. Let 



f --- 



Uv-" - 30f ^ 



r ) dt r 



SOf-" + — 
3r 



and r — (r cos 6',r sin 6',0); then, equation (^2) can be written in polar coordinates as 



r9 + 2r9 = Sfg - elCg, 



(Bl) 



where 



fr = /i cos 9 + f2 sin d, /e = -/i sin 6 + f2 cos 6, 



/Cr = r(/CiiCos26' + /Ci2sin26'), /Ce = r(-/Cii sin 261 + /C12 cos 26*). (B2) 

In determining ICr and ICe, we have utilized the fact that the matrix (ICij) is symmetric and traceless. In what follows, 
we let 

Fr := 6fr - elCr, Fg Sfe - elCg (B3) 

for the sake of simplicity. 

The Delaunay elements are defined by the ellipse tangent to the path at time t. The state of relative motion 
is specified by (r,9,Pr,P9) in polar coordinates; therefore, using the standard formulae for elliptic motion we have 
PrPe = esinv andpg/r — l + ecosw, which uniquely specify the eccentricity e and the true anomaly v of the osculating 
ellipse. Next, the semimajor axis of this ellipse a is obtained from pg = a(l — e^) and the eccentric anomaly it can be 
determined from 

cos u — e . ^ \/l — sin u 

cos V = , sm V = . 

1 — e cos u 1 — e cos u 

The Delaunay elements are then given by 

L = y/a, G = pe, 

i^u — esinu, g = 9 — v. (B4) 
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To obtain the equations of motion in Delaunay elements, we observe that the energy and orbital angular momentum 
of the relative motion per unit mass are 



1 

2^' 



(B5) 



where the orientation of the spatial coordinate system is chosen such that pe > 0. Now 

1 dE ^ 
_ =F.v, 

fji at 



where F ■ v = r + rO. It follows that 



dL 



dt ^/l^. 



FrB sinv + Ft 



a(l-e2) 



(B6) 



wher e w e have replaced f by e sinw/pe, which follows from equation (B5). It is a direct consequence of equations (Bl) 
and (^) that 



dG 



(B7) 



To determine d£/dt and dg/dt, we first note the following useful relationships that are obtained from equations (86) 
and (B7), the definitions of the Delaunay elements, and the Kepler law lu = \/L^, 



da 

dt uj 



a(l-e2) 



de VT 



dt 



Consider first dg /dt — d9/dt — dv /dt, where 



Fr e sin v + Fe 



■ - 7^ / r + a 
Fr sm V + Fg e H cos v 



Differentiating the relationship 



d0 _ y/a(l - e2) 
di ^ ^ 



In r — In a + In (1 — e^) — In (1 + e cos i)) 



(B8) 



(B9) 



with respect to time and using the definition of g together with equations (B8) and ( |B9| ), we have 

dg \/ a(l — e^) 



dt 



'Fr cosv + Fg (l + ——^ — ^ ) sinw 
V a(l — e'^)/ 



(BIO) 



Let us now compute d£/dt using the definition of i: The Kepler equation £ = u — esinw expresses the mean anomaly 
£ in terms of the eccentric anomaly and the eccentricity of the orbit. Thus, 



d£ r du 



dt a dt avT" 



de 



(Bll) 



Using the identity 



r a — r er sin v du 

r — a — he = . — , 

a e VI - at 



which is obtained by differentiating r — a(l — ecosw) with respect to time, and equation (B8), we find 
d£ 



dt eJa 



\Fr (— 2e + cosu + e cos^ v) — Fg {2 ^ e cosv) sin-C] 



(B12) 
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Therefore, the eq uation s of m otion in Cartesian coordinates are equivalent to Delaunay's system given by equa- 
tions ^^), (|b^), (B12), and ( |B1CI|) . These equations for the planar problem are closely related to corresponding 



equations for the variation of the orbital elements of the osculating ellipse. In the general, i.e. three-dimensional, 
Kepler problem, the osculating ellipse is characterized by six orbital elements and the variation of these orbital el- 
ements due to perturbing forces is given by a set of equations that are usually associated in celestial mechanics 
with the name of Lagrange (e.g. "Lagrange's planetary equations"). In general, the path of the osculating ellipse 
in the six-dimensional manifold of orbital elements encounters singularities at e = and e = 1. These singularities, 
corresponding to a circle (G = L) and a straight line (G = 0), respectively, are also evident in the Delaunay equations. 
The right hand side of the set of Delaunay equations should also be represented in Delaunay elements. To accomplish 



this, recall equations (B3) and let Sfo — eA^n, where D denotes any one of the Delaunay variables, be the net 
contribution of radiation reaction and external forces to the rate of variation of Delaunay element D with respect to 
time. We have shown in a previous paper ||^ that 

^, = - + %m ) , (B13) 



where 



and 



(j){t) = iar^^ cos (m), i!{t) = cos [VLt + p), (B14) 



C{L,GJ,g)^^a^e^cos2g 



+a'^ ^^(v4i, cos 2g cos iy£ — sin 2g sin vH.)^ 

u=l 

5 



S{L,G,i,g)=^-a^^sm2g 



-a^ "^^{Ay sin 2(7 cos vl + cos 2g sin (B15) 



Here we have introduced A^, and B^, which are given in terms of the Bessel functions by 

A, = A^-{2ve{l - e^)Jl{iye) - (2 - e')Mve)), 



B,, = — Vl-e^ (eJ^M - 1^(1 ~ e')J,{ve)). 

Let us now compute the contribution of radiation damping terms to the equations of motion. Using equation (|B2| ) 
and the relationships 

.21, 1 2 G^ 

w = , = + J , 

r a a r r 



we find that 



Aesinv ( 10 L'^G'^ 



18G / 20 i2 5X2^2- 
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The substitution of equation ( B16| ) into the Delaunay equations via equation ( p^ ) results in the set of equations 



presented in equation (|l8|), where e = \J\ ~ / L'^; moreover, these relations involve powers of 1/r as well as sinw 
times powers of 1/r. Therefore, to complete the discussion we need Fourier expansions of these expressions in terms of 
cos v(. and sin similar to those given in equation ( B15D . This can be accomplished by extending the classical methods 
of celestial mechanics described, for instance, in the definitive monograph of Watson p9[ |. The final expressions for 
/u, D G {L,G,i,g\, are not used explicitly in this paper; therefore, we do not present them here for the sake of 
brevity. 

Finally, let us note that the Delaunay equations can now be expressed in the form given in equation (|r; 

APPENDIX C: AVERAGE RATE OF DAMPING 

The average contribution of radiation damping to the equations of motion in Delaunay variables is employed in 
this paper in the calculation of the bifurcation function (Section ^ as well as in the description of averaged dynamics 
(Section ^ . In the light of the results of the previous appendix, the quantity fo , where D is any one of the Delaunay 



variables, is given by equation (17) and can be expressed as a Fourier series 

oo 

fniL, G,£,g) — ao + ^^(oi/ cosvi + by sinvl). (CI) 
We are interested in /o, which is the average of Jo with respect to the "fast" angular variable i.e. 

e2T; 



fD:=^ I ^ fD{L,GJ,g)dL 

Jq 



It is clear that fjj = oq , and this is the quantity that is needed in Section ^ 

The contribution of radiation damping to the bifurcation function B is given by AJ^u, where 

Td - / fD{L,G,iut + i,g)dt. (C2) 



Using the change of v aria ble t — Qt/m + i/n, the resonance condition lo — nVL/m, and the periodicity of with 
respect to I, equation ( |C2| ) can be written as = 27mL^/£). It remains to compute Jd- 
Using the relation d£ = r'^dv/ (a^Vl — e^), we can write 

where the angular brackets denote the average over the true anomaly v. The various averages may be evaluated using 

G" /1\ = 1 + + le\ G'» /4\ = 1 + be' + He", 



^2 / sm v\ 1 „4 / sm v\ 1 



2 V 4 



We find that 
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h = 0, /<, - 0, (C5) 
where the last two equations simply follow from the fact that 

< h(cosv) smv >= 

for every contimious function h. 

For the bifurcation function B, only J^l, ^g, and J^g are needed and these are then given by 

^ f 73 , 37 A 

^G = -27mij(8 + 7e2), 

= 0. (C6) 
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FIG. 1. The plots are for system (^) with parameter values e = 10"*, 5/e — 10~^, a = 2, /3 = 2, p = 0, 
and Q, — .6495190528. The top panel shows L = ^/a versus time for the initial conditions {pr,P0,r,O) equal to 
(.38927, .61085, 2.70900, —2.45677). The middle panel shows L versus time and the bottom panel shows G versus L for the 
initial conditions (.3893, .6109, 2.709, -2.4568). Here, L Ki 1.1547 corresponds to (1 : 1) resonance, 1/L^ = Q. 



FIG. 2. The plot of G versus L for system ( p2| ) with initial conditions {pr,pe, r, 9) equal to (.5, 1, 1, 0) and parameter values 
e = 10"^, S/e = 10"^, a = 2, /? = 2, p = 0, and J7 = 1.299038106. The top panel depicts 157000 points, one after each time 
interval 2-k/Q.. The middle panel is a blow up for iterates 110000-130000, the bottom panel for 140000-155000. 
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